/* program to generate theory plots for the ratio :n_NG/n_G
 *
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include <gsl/gsl_integration.h>
#include "massfunction.h"

int main(int argc, char *argv[]) {
  printf("USAGE: mfn q ftNL name <fNL>\n");

  float fNL=0;
  double calA=8.413e-10;
  if (argv[4]=="a"||argv[4]=="cc") {calA=6.72e-10;}
  if (argv[4]=="ccc") {calA=7.6264e-10;}
  if (argv[6]) {fNL = atof(argv[6]);} printf("fNL: %f\n", fNL);
  mfn(calA, atof(argv[1]), atof(argv[2]), atof(argv[3]), argv[4], fNL);
  

}


void mfn(double A, double deltac, double q, double ftNL, char* name, float fNL) {
  char fname[150], string[500], fnamefNL[150];
  FILE *fp;
  FILE *fpfNL;
  double M, ratio, ratiofNL, ratioM3;
  int which;
  // one could also get the chi2 by fitting to the data right here and include chi2/NDF in the filename--the name will be sufficient to identify the appropriate set of files to look for. //
  int ndf,z;
  //zi=0; if (z==0.5) zi=5; if (z==1.) zi=1; if (z==2.) zi=2;
  float chi2ndf[2], chi2M3ndf[2];
  char Opbase[200];
  //sprintf(Opbase, "/home1/02539/sza5154/work/FeederSim/Output1024");
  sprintf(Opbase, "/home/saroj/Research/feedersim/Theory/Output1024");
  for (z=0; z<3; z++) {
    sprintf(fname, "%s/%s/mfn/mbinz%d.dat.ratio.avg", Opbase, name, z);
    fp = fopen(fname, "r");
  
    float col1, col2, col4, chi2=0., chi2M3=0.;
    ndf = 0;
    while (!feof(fp)){
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g %g", &col1, &col2, &col4);
	if (col1>0) {
	  M=col1*9.65e11; // in h^-1 M_sun
	  ratioM3 = ratio3(M, A, deltac, q, ftNL, z);
	  ratio = ratio3(M, A, deltac, q, ftNL, z)+ratio4(M, A, deltac, q, ftNL,z)+ratio5(M,A,deltac,q,ftNL,z);
	  
	  chi2M3+=pow((ratioM3-col2)/col4,2);
	  	  
	  chi2+= pow((ratio-col2)/col4,2);
	  ndf++;
	  printf("r3: %g, r:4 %g, r4f: %g, r5: %g\n", ratio3(M, A, deltac, q, ftNL, z), ratio4(M, A, deltac, q, ftNL,z), ratio4f(M, A, deltac, q, ftNL,z), ratio5(M, A, deltac, q, ftNL,z));

	}
      }
    }
    chi2ndf[z]=chi2/ndf;
    chi2M3ndf[z]=chi2M3/ndf;
    printf("chi2: %g, ndf: %d, chi2/ndf: %g\n", chi2, ndf, chi2ndf[z]);
    printf("chi2M3: %g, ndf: %d, chiM32/ndf: %g\n", chi2M3, ndf, chi2M3ndf[z]);
    fclose(fp);
  /*
    // now write the theory data file along with the chi2/ndf as the last number before the .txt on the filename
    for (which=3; which<5; which++) {
      if (which==3) {
    sprintf(fname, "%s/%s/mfn/theory_%f_z%d_%d_%f.txt", Opbase, name, deltac,z, which, chi2M3/ndf);sprintf(fnamefNL, "%s/%s/mfn/ratio_%f_%f_z%d_%d.txt", Opbase,name, fNL,deltac,z,3);} else {
       sprintf(fname, "%s/%s/mfn/theory_%f_z%d_%d_%f.txt", Opbase, name, deltac,z, which, chi2/ndf);sprintf(fnamefNL, "%s/%s/mfn/ratio_%f_%f_z%d_%d.txt", Opbase, name, fNL,deltac,z,4);}
    fp = fopen(fname, "w");
    
    fpfNL=fopen(fnamefNL,"w");    
    
    for (M=1e12;M<1e16;M=M*1.02) {
      if (which==3) {ratio = ratio3(M, A, deltac, q, ftNL, z);} else {
      ratio = ratio3(M, A, deltac, q, ftNL, z)+ratio4f(M, A, deltac, q, ftNL,z)+ratio5(M,A,deltac,q,ftNL,z);}
      
      fprintf(fp, "%g %g\n", M, ratio);
      if (fNL) {

	if (which==3) {ratiofNL=ratio3(M,AC, deltac, 1, fNL, z);}
	else {ratiofNL=ratio3(M,AC,deltac,1,fNL,z)+ratio4(M,AC,deltac,1,fNL,z);
	}
        fprintf(fpfNL, "%g %g\n", M, ratiofNL);
      }
      
    }
    fclose(fp);
    fclose (fpfNL);
    }
    
    // write data with incremental which from 3 to 6
    for (which=3; which<8; which++) {
      sprintf(fname, "/home/saroj/Research/TwoField/Simulation/Theory/mfns/mfnw_%s_%d_%d", name, z, which);
      fp = fopen(fname, "w");
      for (M=1e12; M<1e16; M=M*1.02) {
	ratio =0;
	if (which>6) ratio+=ratio7(M, A, deltac, q, ftNL, z);
	if (which>5) ratio+=ratio6(M, A, deltac, q, ftNL, z);
	if (which>4) ratio+=ratio5(M, A, deltac, q, ftNL, z);
	if (which>3) ratio+=ratio4(M, A, deltac, q, ftNL, z);
	if (which>2) ratio+=ratio3(M, A, deltac, q, ftNL, z);
	fprintf(fp, "%g %g\n", M, ratio);
      }
      fclose(fp);
    }
    */
  }  
    
    /*
  // once these are done generate gnuplot files that plot simulation data/theory curves 
  sprintf(fname, "%s/%s/mfn/gplot_%f.txt", Opbase, name, deltac);
  fp = fopen(fname, "w");
  fprintf(fp, "set logscale x\n");
  fprintf(fp, "set xrange [2e13:2e15]\n");
  fprintf(fp, "set yrange [0:6]\n"); // can calculate the max (later)
  fprintf(fp, "set size square\n");
  fprintf(fp, "set xlabel \"Mass ($h^{-1}M_{\\\\odot}$)\"\n");
  fprintf(fp, "set ylabel \"$\\\\frac{n_{NG}}{n_G}$\"\n");
  fprintf(fp, "set title \"$\\\\texttt{%s}$\"\n",name);
  fprintf(fp, "plot \"theory_%f_z%d_%d_%f.txt\" title \"$z=%d, \\\\chi^2/n=%.2f$\" with lines lt 1 lc 1, \"mbinz%d.dat.ratio.avg\" using ($1*9.65e11):2:4 notitle with errorbar lt 1 lc 1", deltac,0,4, chi2ndf[0],0,chi2ndf[0],0);
  fprintf(fp, ", \"theory_%f_z%d_%d_%f.txt\" title \"$z=%d, \\\\chi^2/n=%.2f$\" with lines lt 2 lc 3, \"mbinz%d.dat.ratio.avg\" using ($1*9.65e11):2:4 notitle with errorbar lt 2 lc 3", deltac,1, 4, chi2ndf[1],1,chi2ndf[1],1);
  fprintf(fp, ", \"theory_%f_z%d_%d_%f.txt\" notitle with line lt 4 lc 1 lw 0.25, \"theory_%f_z%d_%d_%f.txt\" notitle with line lt 4 lc 3 lw 0.25\n",deltac, 0, 3, chi2M3ndf[0], deltac, 1, 3, chi2M3ndf[1]);
  fprintf(fp, "set term epslatex color\n");
  fprintf(fp, "set output \"%sz01_%f.tex\"\n", name, deltac);
  fprintf(fp, "replot");
  
  fclose (fp);
  */
}

double ratio3(double M, double A, double deltac, double q, double ftNL, double z) {
  return 1 + M3R(M,A,q,ftNL)*He3(nuc(deltac, M, A, q, ftNL,z))/6-M3RD(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/(6*nucD(deltac,M,A,q,ftNL,z));
}

double ratio4f(double M, double A, double deltac, double q, double ftNL, double z) {
  return M4R(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/24-M4RD(M,A,q,ftNL)*He3(nuc(deltac,M,A,q,ftNL,z))/(24*nucD(deltac,M,A,q,ftNL,z));
}

double ratio4(double M, double A, double deltac, double q, double ftNL, double z) {
  double nucp = nucD(deltac,M,A,q,ftNL,z);
  return M2R(M,A,q,ftNL)*He2(nuc(deltac,M,A,q,ftNL,z))/2+M4R(M,A,q,ftNL)*He4(nuc(deltac,M,A,q,ftNL,z))/24+pow(M3R(M,A,q,ftNL),2)*He6(nuc(deltac,M,A,q,ftNL,z))/72-M2RD(M,A,q,ftNL)*He1(nuc(deltac,M,A,q,ftNL,z))/(2*nucp)-M4RD(M,A,q,ftNL)*He3(nuc(deltac,M,A,q,ftNL,z))/(24*nucp)-M3R(M,A,q,ftNL)*M3RD(M,A,q,ftNL)*He5(nuc(deltac,M,A,q,ftNL,z))/(36*nucp);
}

double ratio5(double M, double A, double deltac, double q, double ftNL, double z)
{
  return M5R(M, A, q, ftNL)*He5(nuc(deltac, M, A, q, ftNL, z))/120;
}

double ratio6(double M, double A, double deltac, double q, double ftNL, double z) {
  return (M6R(M, A, q, ftNL)/720 + 0.5*pow(M3R(M, A, q, ftNL)/6,2) )*He6(nuc(deltac, M, A, q, ftNL, z));
}

double ratio7(double M, double A, double deltac, double q, double ftNL, double z) {
  return (M3R(M, A, q, ftNL)*M4R(M, A, q, ftNL)/(48*3)+M7R(M, A, q, ftNL)/(720*7))*He7(nuc(deltac, M, A, q, ftNL, z));
}

double nuc(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {    gf = 0.473345;} else if (z==2) {gf = 0.327534;} else { gf=0.; printf("no growth function value!\n"); }
  return deltac*0.76001/(sqrt(d2R(M,A,q,ftNL))*gf);
}

double nucD(double deltac, double M, double A, double q, double ftNL, double z) {
  double gf;
  if (z==0){gf = 0.76001;} else if (z==0.5) {gf = 0.59455;} else if (z==1) {    gf = 0.473345;} else if (z==2) {gf = 0.327534;} else {  gf=0.; printf("no growth function value!\n"); }
  return -deltac*0.76001*(0.5*(A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
   pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),0.5)/d2R(M,A,q,ftNL)/gf;
}

double M2R(double M, double A, double q, double ftNL) {
  return d2R(M,A,q,ftNL)/d2R(M,AC,0,0)-1;
}

double M3R(double M, double A, double q, double ftNL) {
  return d3R(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),1.5);
}

double M4R(double M, double A, double q, double ftNL) {
  // assume feeder for q <0.1
  if (q<0.1) {
    return 0.3*48 * pow(M3R(M, A, q, ftNL)/8.0, 4.0/3.0);
  }
  return 1.75* d4R(M,A,q,ftNL)/pow(d2R(M,A,q,ftNL),2);
}

double M5R(double M, double A, double q, double ftNL) {
  // assume feeder scaling
  if (q<0.1) {return 2.0* pow(M3R(M, A, q, ftNL)/8.0, 5.0/3.0)*24*16;}
  return An(5)*pow(M3R(M,A,q,ftNL)/8.0,3.0);
}

double An(double n) {return tgamma(n-1.0)*pow(2, n-3);}
double Bn(double n) {return tgamma(n)*pow(2,n-1);}

double M6R(double M, double A, double q, double ftNL) {
  // assume feeder scaling
  return 1*Bn(6.)*pow(M3R(M, A, q,ftNL)/8.0,2.0);
}

double M7R(double M, double A, double q, double ftNL) {
  return 2.5*Bn(7.)*pow(M3R(M, A, q, ftNL)/8.0, 7./3);
}

double M2RD(double M, double A, double q, double ftNL) {
  return -(((A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M))*I2R1D(M))/
      (AC*pow(I2R1(M),2))) + (A*I2R1D(M) + 
      pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M))/(AC*I2R1(M));
}

double M3RD(double M, double A, double q, double ftNL) {
  return A*A*ftNL*q*(-(12*A*(0.75*I3R1(M)+A*ftNL*ftNL*q*I3R2(M))*(I2R1D(M)+A*ftNL*ftNL*I2R2D(M)))/pow(q*q+A*I2R1(M)+A*A*ftNL*ftNL*I2R2(M),2.5)+(6*I3R1D(M)+8*A*ftNL*ftNL*q*I3R2D(M))/pow(q*q+A*I2R1(M)+A*A*ftNL*ftNL*I2R2(M),1.5)); //check this with MM setup.nb numerical results!
}

double M4RD(double M, double A, double q, double ftNL) {
  return (-2*(48*pow(A,3)*pow(ftNL,2)*pow(q,3)*I4R1(M) + 
        48*pow(A,4)*pow(ftNL,4)*pow(q,4)*I4R2(M))*
      (A*I2R1D(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2D(M)))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),3) + 
   (48*pow(A,3)*pow(ftNL,2)*pow(q,3)*I4R1D(M) + 
      48*pow(A,4)*pow(ftNL,4)*pow(q,4)*I4R2D(M))/
    pow(A*I2R1(M) + pow(A,2)*pow(ftNL,2)*pow(q,2)*I2R2(M),2);
}

double d2R(double M, double A, double q, double ftNL) {
  return d2R1(M, A, q, ftNL) + d2R2(M, A, q, ftNL);
}

double d3R(double M, double A, double q, double ftNL) {
  return d3R1(M, A, q, ftNL) + d3R2(M, A, q, ftNL);
}

double d4R(double M, double A, double q, double ftNL) {
  return d4R1(M, A, q, ftNL) + d4R2(M, A, q, ftNL);
}

double d2R1(double M,double A,double q,double ftNL){return A*I2R1(M);}
double d2R2(double M,double A,double q,double ftNL){return pow(q*ftNL*A,2)*I2R2(M);}
double d3R1(double M,double A,double q,double ftNL){return 6*q*q*ftNL*A*A*I3R1(M);}
double d3R2(double M,double A,double q,double ftNL){return 8*pow(q*ftNL*A,3)*I3R2(M);}
double d4R1(double M,double A,double q,double ftNL){return 48*q*A*pow(q*ftNL*A,2)*I4R1(M);}
double d4R2(double M,double A,double q,double ftNL){return 48*pow(q*ftNL*A,4)*I4R2(M);}

double I2R1(double M) {
  return pow(M, -2.83160)*exp(246.051+58591./pow(log(M),2)-6140.93/log(M));
}

double I2R2(double M) {
  return pow(M, -3.13172)*exp(271.811+64250.2/pow(log(M),2)-6746.78/log(M));
}

double I3R1(double M) {
  return pow(M, -4.31865)*exp(373.355+88145.1/pow(log(M),2)-9265.59/log(M));
}

double I3R2(double M) {
  return pow(M, -7.48587)*exp(660.953+175586./pow(log(M),2)-17909./log(M));
}

double I4R1(double M) {
  return  pow(M, -6.92941)*exp(610.192+155357./pow(log(M),2)-15919.7/log(M));
}

double I4R2(double M) {
  return I4R1(M)*I3R2(M)/I3R1(M); // this is an assumption!
}

double I2R1D(double M) {
  return (-2.8315989340868626*exp(246.05052999073197 + 58590.98533491765/pow(log(M),2) - 6140.926322745515/log(M)))/pow(M,3.8315989340868626) + 
   (exp(246.05052999073197 + 58590.98533491765/pow(log(M),2) - 6140.926322745515/log(M))*(-117181.9706698353/(M*pow(log(M),3)) + 6140.926322745515/(M*pow(log(M),2))))/pow(M,2.8315989340868626);
}

double I2R2D(double M) {
  return (-3.1317214228523023*exp(271.81104558816355 + 64250.21701980513/pow(log(M),2) - 6746.78180491083/log(M)))/pow(M,4.131721422852302) + 
   (exp(271.81104558816355 + 64250.21701980513/pow(log(M),2) - 6746.78180491083/log(M))*(-128500.43403961026/(M*pow(log(M),3)) + 6746.78180491083/(M*pow(log(M),2))))/pow(M,3.1317214228523023);
}

double I3R1D(double M) {
  return (-4.318649684650349*exp(373.3549059620789 + 88145.12900844424/pow(log(M),2) - 9265.593589986682/log(M)))/pow(M,5.318649684650349) + 
   (exp(373.3549059620789 + 88145.12900844424/pow(log(M),2) - 9265.593589986682/log(M))*(-176290.2580168885/(M*pow(log(M),3)) + 9265.593589986682/(M*pow(log(M),2))))/pow(M,4.318649684650349);
}

double I3R2D(double M) {
  return (-7.485871213528042*exp(660.9532909886348 + 175585.91029280794/pow(log(M),2) - 17908.95641931269/log(M)))/pow(M,8.485871213528043) + 
   (exp(660.9532909886348 + 175585.91029280794/pow(log(M),2) - 17908.95641931269/log(M))*(-351171.8205856159/(M*pow(log(M),3)) + 17908.95641931269/(M*pow(log(M),2))))/pow(M,7.485871213528042);
}

double I4R1D(double M) {
  return (-6.92941234239259*exp(610.1924526222084 + 155356.68340157234/pow(log(M),2) - 15919.739557836265/log(M)))/pow(M,7.92941234239259) + 
   (exp(610.1924526222084 + 155356.68340157234/pow(log(M),2) - 15919.739557836265/log(M))*(-310713.3668031447/(M*pow(log(M),3)) + 15919.739557836265/(M*pow(log(M),2))))/pow(M,6.92941234239259);
}

double I4R2D(double M) {
  return (-10.096633871270283*exp(897.7908376487643 + 242797.46468593605/pow(log(M),2) - 24563.102387162275/log(M)))/pow(M,11.096633871270283) + 
   (exp(897.7908376487643 + 242797.46468593605/pow(log(M),2) - 24563.102387162275/log(M))*(-485594.9293718721/(M*pow(log(M),3)) + 24563.102387162275/(M*pow(log(M),2))))/pow(M,10.096633871270283);
}
// the Hermite polynomials

double He1(double x) {return x;}
double He2(double x) {return x*x-1;}
double He3(double x) {return pow(x,3)-3*x;}
double He4(double x) {return pow(x,4)-6*x*x+3;}
double He5(double x) {return pow(x,5)-10*pow(x,3)+15*x;}
double He6(double x) {return pow(x,6)-15*pow(x,4)+45*x*x-15;}
double He7(double x) {return pow(x,7)-21*pow(x,5)+105*pow(x,3)-105*x;}


// the growth function
/*
double hubble(double z) {
  return 100*h*sqrt(Omm*pow(1+z,3)+(1-Omm));
}

double growth(double z) {
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
  double result, error;
  //double alpha=1.0;
  gsl_function F;
  F.function=&gr_int;
//F.params=&alpha;
  gsl_integration_qags(&F, z, 1000, 0, 1e-5, 1000, w, &result, &error);
  return (5*Omm*hubble(z)/(2*100*h))*result;
}

double gr_int(double z) {
  return (1+z)/pow(hubble(z)/(100*h),3);
}
*/
